Effective interaction between helical bio-molecules 
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The effective interaction between two parallel strands of helical bio-molecules, such as deoxyri- 
bose nucleic acids (DNA), is calculated using computer simulations of the "primitive" model of 
electrolytes. In particular we study a simple model for B-DNA incorporating explicitly its charge 
pattern as a double-helix structure. The effective force and the effective torque exerted onto the 
molecules depend on the central distance and on the relative orientation. The contributions of non- 
linear screening by monovalent counterions to these forces and torques are analyzed and calculated 
for different salt concentrations. As a result, we find that the sign of the force depends sensitively 
on the relative orientation. For intermolecular distances smaller than 6A it can be both attractive 
and repulsive. Furthermore we report a nonmonotonic behaviour of the effective force for increasing 
salt concentration. Both features cannot be described within linear screening theories. For large 
distances, on the other hand, the results agree with linear screening theories provided the charge of 
the bio-molecules is suitably renormalized. 

PACS: 87.15.Kg, 61.20Ja, 82.70.Dd, 87.10+e 
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I. INTRODUCTION 

Aqueous solutions of helical bio-molecules like deoxyri- 
bose nucleic acids (DNA) are typically highly charged 
such that electrostatic interactions play an important 
role in many aspects of their structure and function 
Understanding the total effective interaction be- 
tween two helical molecules is important since this gov- 
erns the self-assembly of bio-molecules, like bundle for- 
mation and DNA condensation or compaction which in 
turn is fundamental for gene delivery and gene therapy. 
In aqueous solution, such rod-like polyelectrolytes release 
counterions in the solution which ensure global charge 
neutrality of the system. Together with these counteri- 
ons, there are, in general, added salt ions dissolved in the 
solution. The thermal ions screen the bare electrostatic 
interactions between the bio-molecules, such that the ef- 
fective interaction between them is expected to become 
weaker than the direct Coulomb repulsion. For very high 
concentrations of bio-molecules or short distances even 
a mutual attraction due to countcrion " overscreening" is 
conceivable 

In this paper, we study the effective interaction be- 
tween two parallel helical bio-molecules. In particular, 
we investigate how the electrostatic interactions are in- 
fluenced by details of the charge pattern on the biological 
macromolecules. In fact, in many cases, as e.g. for DNA 
molecules, the charge pattern on the molecules is not 
uniform but exhibits an intrinsic helix structure. If two 
parallel helical molecules arc nearby, this helix structure 
will induce an interaction that depends on the relative 
orientation of the two helices. Our studies are based on 
computer simulation of the "primitive" model of elec- 
trolytes (25). In particular we study a simple model for 
B-DNA. This model explicitly takes into account the 



double-helical charge pattern along the DNA-strand, it 
also accounts for the molecular shape by modeling the 
major and minor grooves along the strand. The charged 
counter- and salt ions in the solutions are explicitly in- 
corporated into our model. On the other hand, the wa- 
ter molecules only constitute a continuous background 
with a dielectric constant e screening the Coulomb in- 
teractions. Hence the discrete nature of the solvent is 
neglected as well as more subtle effects as image charges 
induced by dielectric discontinuities at the DNA-water 
boundary p6|-|29|], hydration effects due to the affection 
of the hydrophilic surface to the interfacial layers of wa- 
ter pfj|-p5[|, and spatial dependent dielectric constants 
resulting from the decreasing water mobility in confining 
geometries and from saturation effects induced by water 
polarization near the highly charged molecular surfaces 

mm. 

Our motivation to consider such a simple "primitive" 
model is threefold: First, though solvent effects seem to 
be relevant they should average out on a length scale 
which is larger than the range of the microscopic sizes. 
Hence the electrostatic effects are expected to dominate 
the total effective interactions. Second, it is justified to 
study a simple model completely and then adapt it by 
introducing more degrees of freedom in order to better 
match the experimental situation. Our philosophy is in- 
deed to understand the principles of a simple model first 
and then turn step by step to more complicated models. 
Third, even within the "primitive" approach, there are 
many unsolved problems and unexpected effects such as 
mutual attraction of equally charged particles. Our com- 
puter simulation method has the advantage that "exact" 
results are obtained that reflect directly the nature of 
the model. Hence we get rid of any approximation in- 
herent in a theoretical description. Consequently, the 
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dependence of the effective interactions on a model pa- 
rameter can systematically be studied and the trends can 
be compared to experiments. In this respect our model is 
superior to previous studies that describe the counterion 



screening by linear Debye-Hiickel fl3 
ear Poisson-Boltzmann theory [|2 




or nonlin- 



- 51| and even to 
recent approaches that include approximatively counte- 
rion correlations |5^,Q. We also emphasize that one 
main goal of the paper is to incorporate the molecular 
shape and charge pattern explicitly which is modelled in 
many studies simply as a homogeneously charged cylin- 
der fnililii. In fact we find that the double - helix 

structure has an important influence on the effective in- 
teraction for surface-to-surface separations smaller than 
6A. In detail, the interaction can be both repulsive and 
attractive depending on the relative orientation and the 
mutual distance between two parallel DNA strands. This 
effect which is typically ignored in the charged-cylinder 
model for DNA will significantly affect the self-assembly 
of parallel smectic layers of DNA fragments and may re- 
sult in unusual crystalline structures at high concentra- 
tions. 

Let us also mention that many theoretical studies in- 
volve only a single DNA molecule j5^-j5^^|. To extract 
the effective interaction, however, one has at least to in- 
clude two molecules in the model which is the purpose of 
the present paper. In this study we only consider mono- 
valent counterions. Multivalent counterions and a more 
detailed survey on the influence of model parameters on 
the effective interactions will be considered in a subse- 
quent publication. 

The remainder of paper is organized as follows. In 
chapter II, we present the details of the model used in this 
paper. Chapter III describes the target quantities of the 
applied model. Simulation details are presented in chap- 
ter IV. Theories based on linear screening approaches 
such as the homogeneously charged cylinder model, the 
Yukawa segment model and the Kornyshev-Leikin theory 
fiof are shortly discussed in chapter V. Results of the sim- 
ulation and their comparison to linear screening theories 
are contained in sections VI- VIII for the point-charge 
model, the grooved model and added salt respectively. 
We conclude in section IX. 



II. THE MODEL 

The charge pattern and the shape of a single B-DNA 
molecule is basically governed by the phosphate groups 
which exhibit a double helix structure with right-hand 
hclicity. We model this by an infinitely long neutral hard 
cylinder oriented in z direction with additional charged 
hard spheres whose centers are located on top of the 
cylindrical surface. Each charged sphere describes a 
phosphate group and hence the spheres form a double 
helix structure. In detail, the effective cylindrical diam- 
eter D is commonly chosen to be D = 20 A |6l|,|62|,f49| . 



The spheres are monovalent, i.e. their charge q p < 
corresponds to one elementary charge e > 0, q p = — e, 
and they have an effective diameter d p . We do not fix 
dp but keep it as an additional (formal) parameter in 
the range between d p = 0.2A (practically the point-like 
charge limit) to d p = &A (to incorporate a groove ge- 
ometry for the molecule). Furthermore, the helical pitch 
length is P = 34A; the number of charged spheres per 
pitch length (or per helical turn) is 10. Consequently, 
successive charges on the same strand are displaced by 
an azimuthal angle of 36° corresponding to a charge spac- 
ing of 3AA in z direction. In a plane perpendicular to the 
z direction, phosphate groups of the two different helices 
are separated by an azimuthal angle of <f> s = 144°, see 
Figure T], fixing the minor and the major helical groove 
along the DNA molecule. 
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FIG. 1. A schematic picture explaining the positions of 
DNA molecules and the definition of the different azimuthal 
angles cf>o, 4>, <f> a . For further information see text. 

We place the discrete charges on the two different he- 
lices such that two of them fall in a common plane per- 
pendicular to the z axis, see again Figure [J The to- 
tal line charge density along the DNA molecule is then 
A = -Q.59e/A. 

The second DNA molecule is considered to be paral- 
lel to the first one in our studies. The separation be- 
tween the two cylinder origins is R, we also introduce the 
surface-to-surface separation h = R — D. The position 
of the two double helices can be described by a relative 
angle difference <fi between the two azimuthal angles de- 
scribing the position of the bottom helix with respect to a 
fixed axis in the xy plane. This is illustrated in Figure [|. 
The relative orientation <f> is the key quantity in describ- 
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ing the angle dependence of the forces induced by the 
helical structure. We remark that we only study a situ- 
ation where the discrete phosphates from different DNA 
strands possess the same z coordinates for <j> = 0. Small 
shifts in the z coordinate are not expected to change the 
results significantly. A further parameter characterizing 
the discrete location of the phosphate charges along the 
strands is the azimuthal angle 4>q of a phosphate charge 
with respect to the cylinder separation vector, see again 
Figure All results are periodic in cj>o with a periodicity 
of 36°. 

In addition to the DNA molecules we describe the 
countcrions by charged hard spheres of diameter d c and 
charge q c . The counterions are held at room temper- 
ature T = 298K. Their concentration is fixed by the 
charge of the DNA molecules due to the constraint of 
global charge neutrality. Also, additional salt ions with 
charges q+ and q- , modelled as charged hard spheres of 
diameters d + and d_, are incorporated into our model. 
The salt concentration is denoted by C s . The discrete 
nature of the solvent, however, is neglected completely. 

The interactions between the mobile ions and phos- 
phate charges are described within the framework of the 
primitive model as a combination of excluded volume and 
Coulomb interactions screened by the dielectric constant 
e of the solvent. The corresponding pair interaction po- 
tential between the different charged hard spheres is 



Vn{r) 



DO 



for r < (di +dj)/2 
for r > (di + dj)/2 



(1) 



where r is the interparticle separation and i,j are indices 
denoting the different particles species. Possible values 
for i and j are c (for counterions), +, — (for positively 
and negatively charged salt ions), and p (for phosphate 
groups). In addition, there is an interaction potential 
between the DNA hard cylinder and the free ions 
i = c, +, — which is of simple excluded volume form such 
that these ions cannot penetrate into the cylinder. 

Due to the length of this paper and the large number of 
quantities, we summarize most of our notation in Table |. 



III. TARGET QUANTITIES 

Our target quantities are equilibrium statistical av- 
erages for the local counter- and salt ion densities and 
the effective forces and torques exerted onto the bio- 
molecules. For that purpose we consider a slightly more 
general situation with N parallel DNA molecules con- 
tained in a system of volume V. The cylinder centers are 
fixed at positions B4 (i = 1, ...,7V) in the rry-plane. We 
further assume that there are iV c counterions and N + , A_ 
salt ions in the same system. By this we obtain partial 
concentrations n c — N c /V,n + = N + /V,ri- = N-/V of 
counter and salt ions. 



First we define the equilibrium number density profiles 
Pj(r) (j = c, + , — ) of the mobile ions in the presence of 
the fixed phosphate groups via 



N 3 



(2) 



Here {rf} denote the positions of the ith particle of 
species j. The canonical average < ... > over an {in- 
dependent quantity A is defined via the classical trace 



(A) 



N N + N- 

k=\ J m=l J n=l J 



exp(-/3 ]T [V°+ J2 Uij})xA 



(3) 



=c, + , 



]=c,p, + ,- 



Here (3 = 1/ksT is the inverse thermal energy (fcs de- 
noting Boltzmann's constant) and 



^ = (i-^)EE^i r l-^i) 



(4) 



1=1 k=l 



is the total potential energy of the counter- and salt ions 
provided the phosphate groups are at positions {r^} (n = 
1, ...,N p ). Finally the prefactor 1/Z in eq.(||) ensures 
correct normalization, < 1 >= 1. Note that the density 
profiles pj (r) also depend parametrically on the positions 
{t%} of all the fixed phosphate groups (n = 1, N p ). 

Now we define the total effective force Fi per pitch 
length acting onto the ith DNA molecule (i = 1, ...,N). 
As known from earlier work [|63U64|Jlll , |65[| it contains three 
different parts 



Pi = 



The first term, , is the direct Coulomb force acting 
onto all phosphate groups belonging to one helical turn 
of the ith DNA molecule as exerted from the phosphate 
groups of all the other DNA molecules: 



R 



(2) 



R 



(3) 



(5) 



F, 



(i) 



-E 

k 



N p 

E 

n=l; n 7^ fc 



V PP (\f%-f%\) 



(6) 



where the sum J"\ only runs over 10 phosphates belong- 
ing to one helical turn of the ith DNA molecule. This 
term is a trivial sum of direct interactions. 

The second term F£ involves the electric part of 
the interaction between the phosphate groups and the 
counter- and salt ions. Its statistical definition is 



R 



(2) 



=c,+,- 1=1 



rl I)) (7) 



and describes screening of the bare Coulomb interaction 
(P) by the counter and salt ions. 
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TABLE I. List of key variables 


D 


DNA diameter 


d c 


counterion diameter 


dp 


phosphate diameter 


d+, d- 


salt ion diameters 


P 


helical pitch length 


L 


length of simulation box 


£ 


dielectric constant of DNA and water 


T 


temperature 


N p 


number of phosphates in the simulation box 


N c 


number of counterions in the simulation box 


N s 


number of salt ion pairs in the simulation box 


c s 


salt concentration 




counterion valency 


Qp 


phosphate valency 


q+,q- 


salt ion valencies 


X 


linear charge density of the DNA molecule 


Xb 


Bjerrum length 


Tpc 


coupling parameter between phosphates and counterions 


F 


interaction force per pitch length 


F 


used unit for force , Fo = (473) 2 


M 


torque acting onto the DNA molecules 


R 


interaxial separation between DNA molecules 


h 


surface-to-surface separation between DNA molecules 




relative orientational angle between two DNA molecules 


<t>0 


reference orientational angle for one DNA molecule 




interaction force per pitch length within the homogeneously charged cylinder model 




Debye screening length 


p(YS) 


interaction force per pitch length within the Yukawa segment model 


* 


effective phosphate radius in the Yukawa segment model 


<h 


effective phosphate charge in the Yukawa segment model 


C 


size correction factor in the Yukawa segment model 


p(KL) 


interaction force per pitch length within Kornyshev-Leikin theory 





condensation parameter of counterions 



Finally, the third term F£ describes a contact (or de- 
pletion) force arising from the hard-sphere part in V P i (r) 
and V° (i = c, +,-). It can be expressed as an inte- 
gral over the molecular surface Si associated with the 
excluded volume per one helical turn of the ith DNA 
molecule: 



i? (3) = -k B T 



df 



Si 



(8) 



=c,+,- 



where / is a surface normal vector pointing outwards 
the DNA molecule. This depletion term is usually ne- 
glected in any linear electrostatic treatment but becomes 
actually important for strong Coulomb coupling r pc as 
conveniently defined by |ll]j6qj65|| 



2A, 



(9) 



with the Bjerrum length X B = q^e 2 /ek B T. When T pc 
is much larger than one, the Coulomb interaction dom- 
inates thermal interactions and counterion condensation 



may occur. For DNA molecules this is relevant as 
dp + d c = 4 — 6 A and Xb — 7. 14 A for a monovalent 
counterion in water at room temperature, resulting in a 
coupling parameter r pc larger than one. 

Our final target quantity is the total torque per pitch 
length acting onto the ith DNA molecule. Its component 
Mi along the z-direction (with unit vector e z ) can also 
be decomposed into three parts 



Mi = M. 



(i) 



M, 



(2) 



M, 



(3) 



(10) 



with 



E 



N„ 



n— l;n^k 



(11) 



mP = -r. 



E 



i=c,+,- 1=1 



E^i(K-rj|)) 



(12) 
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and 



Mi 3) = k B Te z 



dfx 



(13) 



IV. COMPUTER SIMULATION 

Our computer simulation was performed within a sim- 
ple set-up which is schematically shown in Figure ^. We 
consider two parallel DNA molecules in a cubic box of 
length L with periodic boundary conditions in all three 
directions. L is chosen to be three times the pitch length 
P such that there are N p = 120 phosphate charges in the 
box. The number of counterions iV c = 120 in the box 
is fixed by charged neutrality while the number of salt 
ions, N s , is governed by its concentration C s . The sep- 
aration vector between the centers of the two molecules 
points along the rr-direction of the simulation box. The 
relative orientation is described according to our notation 
presented in chapter II, see again Figure |l} 

We performed a standard Molecular Dynamic (MD) 
code with velocity Verlet algorithm ^57|. System param- 
eters used in our simulations are listed in Table [J. The 
time step At of the simulation was typically chosen to be 
10~ 2 ^Jm dfn/e 2 , with m denoting the (fictitious) mass 
of the mobile ions, such that the reflection of counteri- 
ons following the collision with the surface of DNA core 
cylinder and phosphates is calculated with high precision. 
For every run the state of the system was checked during 
the simulation time. This was done by monitoring the 
temperature, average velocity, the distribution function 
of velocities and total potential energy of the system. On 
average it took about 10 4 MD steps to get into equilib- 
rium. Then during 5 • 10 4 — 5 • 10 6 time steps, we gathered 
statistics to perform the canonical averages for calculated 
quantities. 

The long-ranged nature of the Coulomb interaction 
was numerically treated via the efficient method pro- 
posed by Lekner |^|. A summary of this method is 
given in Appendix A. In order to save CPU time, the 
Lekner forces between pair particles were tabulated in 
a separate code before entering into the main MD cycle. 
The tabulation on a 510 x 510 x 510 grid with spatial step 
=0.1 A was done in the following manner. The first parti- 
cle was fixed at the origin (0,0,0) while the second charge 
was successively embedded on sites of the generated grid. 
Then the force components acting onto the first charge 
were calculated via the Lekner method. A force data 
file was created which was used as a common input for 
all subsequent MD runs. To decrease error coming from 
a finite grid length, the forces in the simulations were 
calculated using the four-step focusing technique ||| . 




FIG. 2. Schematic view of the set-up: Two cylindrically 
shaped DNA molecules with a distance R at positions Ri and 
R2 are placed parallel to the z-axis inside a cube of length L. 
The large gray spheres are counterions of diameter d c . The 
black spheres of diameter d p , connected by the solid line, are 
phosphate charges on the cylindrical surface of diameter D. 
P is the pitch of DNA. Arrays r v and r c point to positions of 
phosphates and counterions. For sake of clarity, the positions 
of added salt ions are not shown. There are periodic boundary 
conditions in all three directions. 



V. LINEAR SCREENING THEORY 

Linear screening theory can be used to get explicit an- 
alytical expressions for the effective interactions between 
helical bio-molecules. These kind of theories, however, 
should only work for weak Coulomb coupling and thus 
represent a further approximation to the primitive model. 
Depending on the form of the fixed charge pattern char- 
acterizing the biomolecules, one obtains different approx- 
imations. 



A. Homogeneously charged cylinder 

The simplest approach is to crudely describe the 
biomoleculc as a homogeneously charged cylinder. In this 
case, the effective interaction force per pitch length be- 
tween two parallel rods reads Eq] 



F = F {HC) = 



2A 2 PA D Ad(r/A g ) r 
e(D/2) 2 Kf(D/(2X D ))r 



(14) 



5 



Typeset using REVTeX 



TABLE II. Parameters used for the different simulation runs. The Debye screening length \p, as defined by Eqn.(|f5|), and 
the Coulomb coupling T pc are also given. 



Run 


d c (A) 


d p (A) 




C S {M) 


A D (A) 


r 

pc 


A 


l 


0.2 


~ 


~ 


9.6 


12 


B 


2 


2 






9.6 


3.6 


C 


2 


6 






9.6 


1.8 


D 


1 


0.2 


15 


0.025 


8.6 


12 


E 


1 


0.2 


60 


0.1 


6.8 


12 


F 


1 


0.2 


120 


0.2 


5.6 


12 


G 


1 


0.2 


440 


0.73 


3.3 


12 


H 


1 


0.2 


1940 


3.23 


1.7 


12 


I 


2 


2 


120 


0.2 


5.6 


3.6 



Here r is the axis-to-axis separation distance between 
cylinders, Ad is the Debye-Hiickel screening length fixed 
by 



ek B T 



4irj(n c {q c e) 2 + n + (q + e) 2 + n_(<7_e) 2 ) 



(15) 



where the factor 7=1 — V cy i/V is a correction due to the 
fact that the mobile ions cannot penetrate into the cylin- 
dric cores which excludes a total volume V cy i. Further- 
more, K\(x) is a Bessel function of imaginary argument. 
Obviously, the torque is zero for this charge pattern. 



B. Yukawa segment model 

It is straightforward to generalize the traditional 
Debye-Hiickel approach to a general charge pattern re- 
sulting in a Yukawa-segment (YS) model f74[ . 
One phosphate charge interacts with another phosphate 
charge via an effective Yukawa potential [frE] 



U(r) = 



(<7pC) 2 



■ exp(-r/A D ) 



(16) 



Here, £ describes a size correction due to the excluded 
volume of the phosphate groups. This term is assumed to 
be of the traditional Derjaguin-Landau-Verwey-Overbeek 
(DLVO) form 



C = exp(r;A D )/(l + r;AD) 



(17) 



where r* = (d p +d c )/2 is an effective phosphate radius for 
the phosphate counterion interaction. We remark that 
nonlinear screening effects and the excluded volume of 
the cylinder can also be incorporated by replacing the 
bare phosphate charge q p with an effective phosphate 
charge q* p 

Using the same notation as in chapter III, the total 
effective force per pitch length acting onto the ith bio- 
molecule is 



Ft = ft 



(YS) 



E ^ 



f%\) (18) 



within in the Yukawa segment model where the sum 
has the same meaning as in Eqn.(||). Note that the con- 
tact term (||) is typically neglected in linear screening the- 
ory. Furthermore, the effective torque per pitch length 
is 



Mi = Mf S) = -e z ■ J2 f% x I V 



N p 

E 

n— l;n^k 



(19) 

There are also analytical expressions for the equilibrium 
density profiles of the mobile ions involving a linear su- 
perposition of Yukawa orbitals around the phosphate 
charges |77 which, however, we will not discuss further 
in the sequel. 



C. Kornyshev-Leikin theory 

The linear Debye-Hiickel screening theory was recently 
developed further and modified to account for dielectric 
discontinuities and counterion adsorption in the grooves 
of the DNA molecule by Kornyshev and Leikin (KL) 
p0| , [78| -|8l]] . An analytical expression for the effective pair 
potential Vkl {R, <j>) per pitch length between two paral- 
lel rods of separation R with relative orientation </> was 
given for separations larger than R > D + A^> . Here we 
only discuss the leading contribution in the special case 
of no dielectric discontinuity which reads 

Vkl{RA) = ^T ^ (-1) £J77 



f3 n f(K' n (k n D/2)) 2 
(20) 



and corresponds to the interaction of helices whose 
strands form continuously charged helical lines. In 
Eqn.pC 



Pn = 



ngK n (k n D/2)l' n (ngD/2) 
k n K' n (k n D/2)I n (ngD/2)' 



(21) 
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k n = y/l/\% + (ng)*, g=^, (22) 

K n and I n are modified Bcsscl functions of nth order, 
and K n (x) = dK n (x)/dx, I n {x) = dl n (x)/dx. 

We emphasize that the KL-theory does not assume a 
priori the double helical phosphate charge pattern as de- 
fined in chapter II. There are rather more possible charge 
patterns considered including a condensation of countcri- 
ons in the minor and major groove along the phosphate 
strands, and on the cylinder as a whole. This involves 
four phenomenological parameters as a further input for 
the KL theory which makes a direct comparison to the 
simulation data difficult. In fact, for the charge pattern 
given in chapter II, the KL-theory reduces to the Yukawa- 
segment model. 

In detail, the charge pattern is characterized by the 
form factor P n 

P n = (1 - A - h ~ h)BS nfi + 

/i0 + / 2 (-l) n 0- (l-/ 3 0) cos(n0 s ). 

Here 8 n , m is the Kronecker's delta function; 9 is the first 
phenomenological input parameter which describes the 
fraction of counterions that are condensed on the whole 
cylinder. The three numbers /, denote the fractions of 
counterions in the middle of the minor groove (/i), in the 
middle of the major groove (/2), and on the phosphate 
strands (f^) with respect to all condensed counterions. 
We note that the sum in ( p(j| ) rapidly converges, such 
that it can safely be truncated for \n\ > 2. It is straight- 
forward to obtain the effective force and torque per pitch 
length between two molecules from (^0|) by taking gradi- 
ents with respect to R and <p. 

VI. RESULTS FOR POINT-LIKE CHARGES AND 
NO ADDED SALT 



first DNA molecule: along a phosphate strand and along 
the minor and major groove. In order to reduce the sta- 
tistical error we course-grained this density field further 
in a finite volume which is illustrated in Figure ||. 




FIG. 3. A schematic picture to explain the procedure of 
counterion density calculations along one pitch length of a 
DNA molecule. The filled circles connected with solid line 
are phosphate groups. The shaded areas correspond to a path 
along the major groove and along one phosphate strand. The 
considered volume has a height £ and width 5. The neigh- 
bouring DNA molecule is assumed to be on the right hand 
side. 



In what follows, we consider the set-up of two parallel 
bio-molecules with periodic boundary conditions shown 
in Figure 2. We projected F\ onto the vector R, defining 
F = Fi ■ (Ri — R2)/ I R\ — R2 |. Hence a negative sign 
of F implies attraction, and a positive sign repulsion. 
The torque is given for the first DNA molecule, hence 
M = Mi. We start with the case of no added salt. First, 
we assume the counterion and phosphate diameters to be 
small, in order to formally investigate the system with a 
high coupling parameter T pc > 10. 

A. Distribution of the counterions around the DNA 
molecules 



This volume is winding around the molecules with a 
height £ and width S. We choose £ = ZAA and 5 — 
2\+d c /2. In Figure [| we plot this coarse-grained density 
field p c (<p) versus the azimuthal angle angle (p from 0° to 
360° where tp is 0° resp. 360° in the inner region between 
the DNA molecules. 

Obviously, the counterion density profile has maxima 
in the neighbourhood of the fixed phosphate charges. 
Furthermore the concentration of counterions is higher 
in the minor than in the major grooves with the in- 
dependence reflecting again the position of the phosphate 
charges. Also in the inner region between the two DNA 
molecules, there are on average more counterions than in 
the outside region. 



We calculated the equilibrium density field (||) of the 
counterions in the vicinity of the DNA molecules by com- 
puter simulation. In detail, we considered three different 
paths to show the counterion density profile around the 
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FIG. 4. Equilibrium counterion density profile p c (<p) in 
units of 1/hDS versus azimuthal angle ip for the parameters 
of run A, <f> = 0° and a rod separation of R — 30A. Solid line: 
counterion density profile along a phosphate strand (due to 
symmetry, the counterion density profiles on the two phos- 
phate strands are the same). Dashed line: counterion density 
profile along the major groove. Dot-dashed line: counterion 
density profile along the minor groove. 



B. Nearly touching configurations 

Let us now consider very small surface-to-surface sep- 
arations between the DNA molecules. In this case one 
expects that the dependence of the forces and torques 
on the relative orientation <f> is most pronounced. For 
such nearly touching configurations, however, the dis- 
creteness of the phosphate charges, as embodied in the 
parameter <po, strongly influences the results as well. The 
qualitative behaviour of the (f> dependence can be under- 
stood from Figure |^. Here two touching DNA molecules 
are shown for different relative orientations <fi where the 
phosphate strands are schematically drawn as continuous 
lines. For certain angles 4> which we call touching angles, 
two neighbouring phosphate charges hit each other. Pos- 
sible touching angles are <j> = 36°, 180°, 324°. If </> is 
chosen to be zero, then two point charges are opposing 
eachother directly. Hence a strong dependence on cj) and 
on <j>Q is expected near touching angles. 

Results from computer simulation and YS-theory are 
presented in Figure |6| The parameters are from run A 
(see Table II) but with d c = 0.8A. The surface-to-surface 
separation is h = 2A. 



36 108 180 252 324 

<)) [degrees] 

FIG. 5. Schematic picture of a DNA-DNA configuration 
for close separation distances. The abscissa corresponds to 
the rotation angle of the first DNA molecule. The second 
DNA molecule is fixed. 

For touching angles, the interaction force becomes 
strongly repulsive. The strongest repulsion is achieved 
for tj> — 180° since two phosphate strands are meeting 
simultaneously. For relative orientations different from 
a touching angle, the force becomes smaller and can be 
both, attractive and repulsive. YS-theory always predicts 
a repulsive force. Again there are strong peaks for touch- 
ing angles in qualitative agreement with the simulation. 
The actual numbers predicted by YS-theory, however, are 
much too large and off by a factor of 6-7 around touching 
angles. 

The torque shows an even richer structure as a func- 
tion of 4>. Near a touching angle it exhibits three zeroes 
corresponding to an unstable minimum exactly at the 
touching angle and two stable minima near the touching 
angles. The YS-theory shows 2 times larger values for 
the torque as compared to the simulation data. 

A qualitatively different force-angle behavior is ob- 
served for a larger counterion diameter. Results for 
d c = lA are shown in Figure 

Here at touching angles, the interaction force is attrac- 
tive. The physical reason for that are the contact forces 
as given by Eqn.(||). Caused by the larger counterion 
diameter, counterions are stronger depleted in the zone 
between the DNA molecules. The torque has qualita- 
tively the same behaviour as before. 

We emphasize that the results do also depend strongly 
on <fio. For 0o = 18°, for instance, the force F practically 
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vanishes for any relative orientation i 
same data for 00 = 0°. 



■ as compared to the 



C. Distance-resolved forces 

We now discuss in more detail the distance-resolved 
effective forces. For the parameters of run A, simulation 
results for F are presented in Figure |^. 

For (f>o = 0, the force depends on the relative orien- 
tation up to a surface-to-surface separation h m 6A 
in accordance with Figure 0. On the other hand, for 
4>o = 18°, there is no dependence at all for any sepa- 
ration. This supports the conclusion of previous works 
|f37| , [55f , that the effect of discreteness of the DNA phos- 
phate charges on the counterion concentration profile is 
small in general and dwindles a few Angstroms from the 
DNA surface. In fact, for h > 6A, there is neither a (f> 
nor a cf> dependence of the force, and the total force is 
repulsive. 

Furthermore we compare our simulation results with 
the prediction of linear screening theories in Figure |^. 
First of all, our simulation data for the total force (solid 
circles) are decomposed into the electrostatic part 
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FIG. 7. Same as Figure o but now for d c — lA. 
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FIG. 6. Interaction force F( left y-axis) and torque M 
(right y-axis) for fixed surface distance h = 2A versus relative 
orientation (f> in degrees. The unit of the force is Fq = (jf,) 2 - 
The solid (dashed) line is the simulation result for F (M) 
while the dot-dashed (dotted) line are data from YS-theory 
for F (M). <f>o is chosen to be zero. The counterion diameter 
is d c = 0.8 A. 



FIG. 8. Effective interaction force F acting onto a DNA 
pair versus the center-to-center distance R. The solid line is 
for <j>o = 18° . In this case there is no significant (^-dependence. 
The meaning of the symbols, that correspond to <j>o = 0, is : 
circles- <j> = 180° , squares- <j> = 36° , triangles- <f> = 45° . 

pi 1 ) .j. p( 2 ) (diamonds) and the contact (or depletion) 
part F^ (open circles). While the latter is strongly re- 
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pulsive, the electrostatic part is attractive such that the 
net force is repulsive. Linear screening theories aim to 
describe the pure electrostatic force only. 

Results for linear screening theories on different lev- 
els are also collected in Figure ^. If one compares with 
the total force, the prediction obtained by a homoge- 
neously charged cylinder is repulsive and off by a factor of 
roughly 1.5. A simulation with a homogeneously charged 
rod yields perfect agreement with linear screening theory 
since the Coulomb coupling is strongly reduced as the rod 
charges are now in the inner part of the cylinder. The 
Yukawa-segment theory is repulsive and off by a factor 
of 3. It is understandable that the YS model leads to 
a stronger repulsion than the charged cylinder model as 
the separation of the phosphate charges in the inner re- 
gion between the DNA molecules is shorter than the rod 
center separation. 



30 



the fraction of counterions which are condensed onto the 
DNA within this shell. The actual value for 5 is some- 
what arbitrary, we first took a microscopic shell of width 
6 = 2.5A as well as S = Xb = 7.1 A. Data for 9 versus the 
rod separation are included in Figure [w] for three differ- 
ent combinations of counterion and phosphate diameters. 
It becomes evident that the fraction 9 of condensed coun- 
terions decreases with the rod distance but saturates at 
large separations. 9 also depends on the size of the coun- 
terions and phosphate charges. If the width of the shell 
S is enhanced towards 5 = Xb = 7.1 A, 9 increases again. 
On the other hand, 9 is independent of the relative orien- 
tation (f>. The actual data are consistent with Manning's 
condensation parameter 82 8^] do = X/\q c \Xs — 0.71 
particularly if the width 6 is taken as one Bjerrum length. 
Our data are also in semiquantitative accordance with 
other computer simulations [ [38f and nuclear magnetic 
resonance (NMR) experiments which show that the con- 
densed counterion fractions are in the range of 0.65 to 
0.85 161 or 0.53 to 0.57 |8lJ45||. 





25 30 35 o 

DNA-DNA distance R[A] 

FIG. 9. Theoretical and simulation results for interaction 
force F versus separation distance R. The unit of the force is 
Fo = (jfj) • The parameters are from run A and <j>o = 18°. 
Symbols: • - simulation data for all DNA rotation angles, 
o - the entropic part F^ 3 \ o - the pure electrostatic part 
(F (1) +F (2) ). Solid line: YS theory. Dot-dashed line: homo- 
geneously charged cylinder model. Dashed line: the predic- 
tions of KL theory with fx = 0.1, f 2 = 0.1, f 3 = 0.7, 6 = 0.71. 

The Kornyshev-Leikin theory requires four counterion 
condensation fractions 9, fx, as an input. We 

have tried to determine these parameters from our sim- 
ulation in order to get a direct comparison without any 
fitting procedure. In order to do so, we introduce a small 
shell around the cylinder of width S and determine 9 as 
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FIG. 10. The condensation parameter 9 versus separation 
distance R. From top to bottom: solid line- run A (d c — lA, 
dp = 0.2A), dot-dashed line-run B (d c = 2 A, d p = 2 A), 
dashed line- run C (d c = 2 A, d p — 6 A). The horizontal line 
at 9 — 0.71 indicates the saturation value at large distances 
for a larger 5 = Ib = 7.1 A. This saturation value is the same 
for run A,B, and C. 

According to our results for the counterion density 
distribution (see Figure ^) we fix the minor and major 
groove fractions to fx = 0.1, fa = 0.1, and the strand 
fraction to / 3 = 0.7. Thus, (1 - fx - f 2 - f 3 ) = 0.1 is the 
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fraction of the condensed counterions which is distributed 
neither on the phosphates strands nor on the minor and 
major grooves. The force in KL theory depends sensi- 
tively on 9 but is rather insensitive with respect to /i, 
/2, fs, and <j). If the Bjerrum length is taken as a width 
for the condensed counterions, 9 — 0.71, then the KL the- 
ory underestimates the total force. If, on the other hand, 
a reduced value of 9 = 0.545 is heuristically assumed, 
then the KL theory reproduces the total force quite well. 

A serious problem of the comparison with linear screen- 
ing theories is that the contact term is not incorporated 
in any theory apart from recent modifications ]8^j64| . In 
fact, one should better compare the pure electrostatic 
part which is attractive in the simulation. Consequently, 
none of the linear screening theories is capable to describe 
the force well. This is due to the neglection of correla- 
tions and fluctuations in linear screening theories. From 
a more pragmatic point of view, however, one may state 
that a suitable charge renormalization leads to quanti- 
tative agreement with the total force. In fact, all three 
theories yield perfect agreement if the phosphate charges 
resp. the condensation parameter 9 is taken as a fit pa- 
rameter. For instance, the YS-model yields perfect agree- 
ment with the simulation for distances larger than 26A if 
in Eqn.(fl6|) a renormalized phosphate charge q* = — 0.6e 
is taken replacing the bare charge q p . But this is still 
unsatisfactory from a more principal point a view. 



VII. RESULTS FOR THE GROOVED MODEL 

The groove structure of DNA is expected to be of in- 
creasing significance as one approaches its surface [ p7| . 
We incorporate this in our model by increasing the phos- 
phate diameter towards d p = 2A (run B) and d p = 6A 
(run C). Results for the condensation parameter 9 are 
shown in Figure [l0|. 9 is decreasing with increasing d p 
since the coupling parameter T pc is decreasing which 
weakens counterion binding to the phosphate groups. 
Also the qualitative shape of the counterion density pro- 
files depends sensitively on the groove nature as can be 
deduced from Figure |ll| as compared to Figure || The 
counterion density along the phosphate strands now ex- 
hibits minima at the phosphate charge positions while it 
was maximal there in Figure |[ Furthermore, the coun- 
terion density in the minor grooves is now higher than 
along the strands due to the geometrical constraints for 
the counterion positions which is similar to results of Ref. 
pf. In fact, recent X-ray diffraction &@ and NMR 



spectroscopy 
mechanics [B3 




)l|,|92fl experiments, as well as molecular 
M[ and Monte Carlo simulations [El sug- 
gest that monovalent cations selectively partition into the 
minor groove. This effect is present also in our simple 
model and can thus already be understood from electro- 
statics and thermostatics. 
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FIG. 11. Same as Figure |i| but now for run C and <p = 45° , 
S = 3A. 

An increasing phosphate and counterion size increases 
the effective forces which is shown in Figure Here, as 
4>q was chosen to be 18°, there is no notable dependence 
on the relative orientation (f>. A similar behavior was ob- 
served in a hexagonally ordered DNA system via Monte 
Carlo calculations |pij| . This is understandable as coun- 
terion screening is becoming less effective. We have tried 
to fit the simulation data using a renormalized charge in 
the YS theory. A good fit was obtained for large sepa- 
rations while there are increasing deviations at shorter 
distances. This is different from our results for small ion 
sizes also shown in Figure ^ where the fit was valid over 
the whole range of separations. The adjustable parame- 
ter q* is shown versus the effective phosphate radius r* 
of the YS model in the inset of Figure [l^. It is increasing 
with increasing r* in qualitative agreement with charge 
renormalization models | p5[ . 

We also note that the physical nature of the electro- 
static part of the interaction force undergoes a trans- 
formation upon decreasing the coupling parameter T pc . 
For strong coupling, r pc — 12 (run A), the electrostatic 
part + F^> is attractive (see Figure ||). For mod- 
erate coupling, r pc = 3.6 (run B), it is nearly zero for 
all distances. Finally, for weak coupling, r pc = 1.8 (run 
C) the electrostatic part is elsewhere repulsive. The en- 
tropic part for these three runs is always repulsive 
and does not undergo a significant change. 
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FIG. 12. Interaction force F versus separation distance R. 
The open circles are simulation data for all relative orien- 
tations <f> with 0o = 18°. From bottom to top: d c — lA, 
dp = 0.2 A (run A); d c = 2 A, d p = 2 A (run B); d c = 2 A, 
d p = 6 A (run C). 

The dashed lines are fits by the YS model. From bottom to 
top: fit for the parameters of run A with q p — — 0.6e; fit 
for the parameters of run B with q p = — 0.75e; fit for the 
parameters of run C with q* — — 0.85e. The inset is the vari- 
ation of the renormalized phosphate charge q p versus effective 
phosphate radius r p * . 



VIII. RESULTS FOR ADDED SALT 

Interactions involving nucleic acids are strongly depen- 
dent on salt concentration. Indeed, the strength of bind- 
ing constants can change by orders of magnitude with 
only small changes in ionic strength . Our simula- 

tions show a similar strong salt impact on the interaction 
force. 

When salt ions are added, there is a competition be- 
tween two effects. The first one is the increasing of the 
direct repulsion between molecules as a consequence of 
delocalizing the adsorbed counterions. The second stems 
from the osmotic pressure of added salt that pushes the 
salt ions to occupy the inner molecular region and to 
screen the DNA-DNA repulsion. As we shall show be- 
low, these two effects result in a novel non-monotonic 
behaviour of the force as a function of salt concentra- 
tion. 
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FIG. 13. Interaction force F acting onto a DNA pair ver- 
sus distance for <f> — 0° and <f>o = 18°. The unit of the force 
is Fa — (^fj) 2 - The solid lines are for increasing salt con- 
centration: 1- run D, 2 - run E, 3 - run F, 4 - run G, 5 - 
run H. Dashed line: reference data without salt from run A. 
The inset shows the force versus salt concentration at fixed 
separation R — 26A. 

Simulation results for F versus distance for increas- 
ing salt concentration are presented in Figure [l^. In our 
simulations, counter and equally charged salt ions are in- 
distinguishable. We take d + = d- = d c , \q + \ = = e. 



It can be concluded from Figure 13 that even a small 
amount of salt ions (line 1, run D, C s — 0.025A/ ) signifi- 
cantly enhances the DNA-DNA repulsion (compare with 
the dashed line corresponding to run A, C s = OM). Upon 
increasing the salt concentration, at large separations, 
h > 1(M, the screening is increased in accordance with 
the linear theory. However, at intermediate and nearly 
touching separations, a non-monotonic behaviour as a 
function of salt concentration is observed as illustrated in 
the inset of Figure [l3l In the inset, the maximum of F oc- 
curs for C s = 0.2M . The physical reason for that is that 
added salt ions first delocalize bound counterions which 
leads to a stronger repulsion. Upon further increasing 
the salt concentration, the electrostatic screening is en- 
hanced again and the force gets less repulsive. In order 
to support this picture we show typical microion config- 
urations and investigate also the fraction of condensed 
counterions as a function of salt concentration. 
Simulation snapshots are given in Figure 



14 



where the 

positions of the mobile ions are projected onto the xy- 
planc. 
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FIG. 14. Two-dimensional microion snapshots projected to 
a plane perpendicular to the helices for <j> = 0°, (j> — 18°, 
R = 30 A. The filled circles are the positions of the counterions 
and positive salt ions, the open circles are the positions for 
the negative salt ions (coions). a - run A, b - run F, c -run G. 

A comparison of the salt-free case (Figure [lija) with 
that of moderate salt concentration (C s = 0.2M, Fig- 
ure [l4]b) reveals that the total number of adsorbed coun- 
terions decreases with increasing C s . Furthermore, for 
C s = 0.2M (Figure |l4]b), there are no coions in the inner 
DNA-DNA region. Thus salt ions do not participate in 
screening. Consequently, the DNA-DNA interaction, due 
to derealization of counterions, will be enhanced. Con- 
trary to that, for C s = 0.73M, (Figure |l|c) the salt co- 
and counterions enter into the inner DNA-DNA region 
and effectively screen the interaction force. 

Further information is gained from the fraction 8 of 
condensed counterions which is plotted as a function of R 
for different salt concentrations C s in Figure 15, We de- 
fine 9 as the ratio of condensed counterions coming from 
the molecules with respect to the total number of counte- 
rions stemming from the molecules. As C s increases, the 
saturation of 9 occurs at smaller distances. In the inset of 
Figure |l5| a non-monotonic behaviour of 9 as a function 
of the added salt concentration is visible which again is 
a clear signature of the scenario discussed above. The 
increase of 9 above a certain threshold of salt concentra- 
tion is mainly due to a counterion accumulation outside 
the grooves. A similar trend was predicted by Poisson- 
Boltzmann |9£| and Monte Carlo J61p7| ] calculations in 
different models. 
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rations. Consequently, the details of the charge pattern 
do not matter for large salt concentrations. 
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FIG. 15. Same as Figure p"o[ but now with added salt. Sym- 
bols: A - run A, • - run D, o - run E, o - run F, * - run G, x 
- run H. The inset shows 9 for fixed distance as a function of 
salt concentration: solid line- for R — 26A; dashed line- for 
R = 30A. 

More details of the forces and the comparison to linear 
screening theories are shown in Figures 16, [l7] and [l^. 
For run F, the different parts of the total force are pre- 
sented in Figure 

m As compared to the salt-free case 
(Figure ^ the pure electrostatic part is again attractive 
but much smaller, while the depletion part is repulsive 
and dominates the total force. All three linear models, 
homogeneously charged cylinder model, YS, and KL the- 
ory, underestimate the force. Note that the KL-theory 
with a 9 parameter corresponding to a width S of one 
Bjerrum length and the homogeneously charged cylinder 
model give the same results. Again with a suitable scal- 
ing of the prefactor by introducing a renormalized phos- 
phate charge q* resp. by fitting the condensed fraction 9, 
one can achieve good agreement with the simulation data 
for distances larger than 24A. The fitting parameter q* 
used for the YS-model is —Lie, while the optimal con- 
densed fraction 9 for the KL-theory is 0.2. The optimal 
renormalized phosphate charge q* is shown versus salt 
concentration in Figure |l7|. Note that the usual DLVO 
size correction factor £ is already incorporated in the in- 
teraction, so what one sees are actual deviations from 
DLVO theory. The renormalized charge q* increases with 
increasing C s which is consistent with the works of Del- 
row et al and Stigter |27| . If one simulates the force 
within the homogeneously charged rod model, one finds 
good agreement with our simulation data for large sepa- 
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FIG. 16. Same as Figure ^ but now for run F and 
4> = 0°, 4>o — 18°. The KL theory was adjusted to 
fi = 0.1, / a = 0.1, h = 0.7,61 = 0.71. The results for KL 
theory and homogeneously charged cylinder models coincide 
exactly. 

We also note that our simulations give no notable de- 
pendence of the force on the relative orientation 4> for 
h > 6A. Only for small separations, h < 6A there is a 
slight dependence in agreement with Ref. ]57| . 

Finally we show the influence of the ion and phosphate 
size on the effective force (for the parameters of run I) 
in Figure The electrostatic part of the force is now 
repulsive but the total force is still dominated by the de- 
pletion part. As far as the comparison to linear screening 
theories is concerned, one may draw similar conclusions 
as for Figure |l6|. The fitting parameter q* needed to de- 
scribe the long-distance behaviour within the YS model 
does not depend sensitively on the phosphate and ion 
sizes. With a suitable scaling of the prefactor one can 
achieve good agreement with the simulation data for dis- 
tance larger than 26A. The fitting parameter q* used 
for the YS-model is —Lie, while the optimal condensed 
fraction 9 for the KL-theory is 0.19. Here again, simu- 
lations of the homogeneously charged cylinder model are 
in good agreement with our results obtained for a double 
stranded DNA molecule. 
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IX. COMMENTS AND CONCLUSIONS 

In conclusion, we have calculated the interaction be- 
tween two parallel B-DNA molecules within a "primi- 
tive" model. In particular, we focussed on the distance- 
and orientation-resolved effective forces and torques as a 
function of salt concentration. Our main conclusions are 
as follows: 
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FIG. 17. Fitted renormalized phosphate charge qp 1 in the 
YS model, versus Debye screening length A_d for runs D-H. 
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FIG. 18. Same as Figure H but now for run I and 
4> = 0°, 4>o = 18°. The KL theory was adjusted to 
fi = 0.1, /a = 0.1, / 3 = 0.7,0 = 0.71. Note that the KL and 
homogeneously charged cylinder models produce the same 



First, the interaction force for larger separations is re- 
pulsive and dominated by microion depiction. The ori- 
entational dependence induced by the internal helical 
charge pattern is short ranged decaying within a typi- 
cal surface-to-surface separation of 6A. For shorter sep- 
arations there is a significant dependence on the rela- 
tive orientation <f> and on the discreteness of the charge 
distribution along the strands. As a function of <fi, the 
force can be both attractive and repulsive. This may 
lead to unusual phase behaviour in smectic layers of par- 
allel DNA molecules. Details of the molecular shape and 
countcrion size are important for small separations as 
well. The torque is relatively small except for small sep- 
arations where it exhibits a complicated (^-dependence. 

Second, as a function of added salt concentration we 
predict a non-monotonic behaviour of the force induced 
by a competition between derealization of condensed 
countcrions and enhanced electrostatic screening. This 
effect can in principle be verified in experiments. 

Third, linear screening theories describe the simula- 
tion data qualitatively but not quantitatively. Having 
in mind that the total force is dominated by the deple- 
tion term which is typically neglected in linear screen- 
ing theory, such theories need improvement. On the 
other hand, the different theories predict the correct long- 
distance behaviour, if a phcnomcnological fit parameter - 
as the renormalized phosphate charge q* for the Yukawa- 
segment model or the condensation fraction 8 for the 
Kornyshev-Leikin model - is introduced. The Yukawa- 
segment model can even predict the orientational depen- 
dence of the force and the torque at smaller distances in 
the case of small counterion and phosphate sizes. Hence, 
a phcnomcnological Yukawa segment model can be used 
in a statistical description of the phase behaviour of many 
parallel DNA strands in a smectic layer. 

Future work should focus on an analysis for divalent 
counterions which are expected to lead to a qualitatively 
different behaviour since the Coulomb coupling is en- 
hanced strongly in this case. Also, one should step by 
step increase the complexity of the model in order to 
take effects such as dielectric discontinuities [38 4l|, |27| , |99|| , 
chemical bindings of counterions in the grooves and dis- 
crete polarizable solvents into account. 
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APPENDIX A: LEKNER SUMMATION METHOD FOR FORCES 



In our simulations we account for the long-range nature of the Coulomb interactions via the effici ent method proposed 
by Lekner |]68[ . This method has been successfully applied to partially periodic systems 14 IOC]. For an assembly of 

— "(c) 

N ions in a central cubic cell of dimension L, the Coulomb force 
repetitions of particle j in the periodic system, is 



exerted onto particle i by particle j, and by all 



?(c) _ HQj 



E 



all cells 1 J 



(Al) 



Because of x, y, z symmetry it is sufficient to consider only one component of the force. For the x-component of the 
force we have 
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Here, Ax = x% — Xj, Ay = yi — yj, Az = z\ — Zj, and Kq(z) is the modified Bessel function o f ze ro order. 
For a pair of particles not aligned parallel to the x-axis, the convergence of the sum in (A2) is fast. Thus an 
evaluation of just 20 terms in the sum is enough to get a part-per-million accuracy. The convergence becomes worse 
when simultaneously |Ay| < S and |Az| < S (5 <C L) for the case m = = n. The number of terms needed in the 
sum for a desired accuracy increases rapidly with increasing S. 

If the particles are aligned parallel to the x-axis such that |Ay| + \Az\ = 0, the sum in ( A2) diverges with m = = n. 
For this particular case Fi X is 
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